Soil Color Contrast

1 Introduction

1.1 Soil Color Contrast Class

1.2 Delta-E

2 Applications

2.1 Setup

library(aqp)
library(soilDB)
library(farver)
library(cluster)
library(ape)
library(colorspace)
library(latticeExtra)
library(knitr)
library(MASS)
library(vegan)

2.2 Core Functions

Soil Survey Technical Note 2

2.2.1 Color Contrast Metrics

Compute several indices of color contrast from 2 vectors of Munsell colors, m1 and m2. Comparisons fully vectorized, e.g. comparisons performed element-wise \(m1_i ~\leftrightarrow~ m2_i\) \(...\) \(m1_n ~\leftrightarrow~ m2_n\).

• \(dH\): change in Munsell hue, see Hue Position section below
• \(dV\): change in Munsell value (absolute value of difference)
• \(dC\): change in Munsell chroma (absolute value of difference)
• \(\Delta{E_{00}}\): CIE2000 Delta-E #1, #2
• \(cc\): soil color contrast class (Soil Survey Technical Note 2)

Results as a data.frame.

# examples
m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', '7.5YR 3/4')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3')

# result is a data.frame
d <- colorContrast(m1, m2)

# check
kable(d, row.names=FALSE)

m1	m2	dH	dV	dC	dE00	cc
10YR 6/3	5YR 3/4	2	3	1	31.112562	Prominent
7.5YR 3/3	7.5YR 4/4	0	1	1	9.522489	Faint
10YR 2/2	2.5YR 2/2	3	0	0	6.959582	Faint
7.5YR 3/4	7.5YR 6/3	0	3	1	30.054654	Distinct

Results summarized in a figure.

# more examples
m1 <- c('10YR 6/3', '7.5YR 3/3', '10YR 2/2', '7.5YR 3/4', '2.5Y 6/8', '5B 4/6')
m2 <- c('5YR 3/4', '7.5YR 4/4', '2.5YR 2/2', '7.5YR 6/3', '10YR 2/1', '5GY 3/4')

# graphical comparison
colorContrastPlot(m1, m2)

Simulated redoximorphic feature colors, contrast classes and \(\Delta{E_{00}}\).

m1 <- paste0('7.5YR 4/', 2:8)
m2 <- rep('10YR 5/2', times=length(m1))
colorContrastPlot(m1, m2, labels = c('F3M', 'MAT'), d.cex = 0.8, col.cex = 0.8)

m1 <- paste0('5Y 4/', 5:1)
m2 <- rep('10YR 3/5', times=length(m1))
colorContrastPlot(m1, m2, labels = c('RMX/FED', 'MAT'), d.cex = 0.8, col.cex = 0.8)

m1 <- c('2.5Y 3/2', '5Y 3/2', '2.5GY 3/2', '5GY 3/2', '7.5BG 3/2', '2.5B 3/2')
m2 <- rep('10YR 3/5', times=length(m1))
colorContrastPlot(m1, m2, labels = c('RMX', 'MAT'), d.cex = 0.8, col.cex = 0.8)

2.2.2 Color Contrast Class

Test rules and cases outlined in Soil Survey Technical Note 2.

10YR 6/3 vs 5YR 3/4

contrastClass(v1=6, c1=3, v2=3, c2=4, dH=2, dV=3, dC=1, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma       res
## 1  6  3  3  4  2  3  1   FALSE   FALSE   FALSE            FALSE Prominent
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3       res
## 1  6  3  3  4  2  3  1   FALSE   FALSE   FALSE Prominent

7.5YR 3/3 vs 7.5YR 4/4

contrastClass(v1=3, c1=3, v2=4, c2=4, dH=0, dV=1, dC=1, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma   res
## 1  3  3  4  4  0  1  1    TRUE   FALSE   FALSE            FALSE Faint
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3   res
## 1  3  3  4  4  0  1  1   FALSE   FALSE   FALSE Faint

10YR 2/2 vs 2.5YR 2/2

contrastClass(v1=2, c1=2, v2=2, c2=2, dH=0, dV=0, dC=0, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma   res
## 1  2  2  2  2  0  0  0    TRUE   FALSE   FALSE             TRUE Faint
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3   res
## 1  2  2  2  2  0  0  0   FALSE   FALSE   FALSE Faint

7.5YR 3/4 vs 7.5YR 6/3

contrastClass(v1=3, c1=4, v2=5, c2=3, dH=0, dV=3, dC=1, verbose = TRUE)

## $faint
##   v1 c1 v2 c2 dH dV dC f.case1 f.case2 f.case3 low.value.chroma      res
## 1  3  4  5  3  0  3  1   FALSE   FALSE   FALSE            FALSE Distinct
## 
## $distinct
##   v1 c1 v2 c2 dH dV dC d.case1 d.case2 d.case3      res
## 1  3  4  5  3  0  3  1    TRUE   FALSE   FALSE Distinct

2.2.3 Hue Position

Differences in Munsell hue (\(dH\)) are computed according to radial position (clock-wise) from 5R    to    5PB.

Use huePosition() to extract hues in order, or convert a vector of hues into an index of positions.

# hues used in the description of soil color (except for neutral hues)
hues <- huePosition(x=NULL, returnHues = TRUE)
# convert hue names into position
hue.pos <- huePosition(hues)
# check
kable(head(cbind(hue.pos, hues)))

hue.pos	hues
1	5R
2	7.5R
3	10R
4	2.5YR
5	5YR
6	7.5YR

Arranged in CIELAB colorspace, the ordering of hues looks like this. The huePositionPlot function in {sharpshootR} will make this kind of figure for any combination of value and chroma.

# make these into real colors by fixing value and chroma at '6'
colors <- paste0(hues, ' 6/6')
# convenient labels for figure
hue.labels <- sprintf("%s\n%s", hues, hue.pos)

# combine colors names, hex representation of colors, and CIELAB coordinates
x <- data.frame(
  colors, 
  hex=parseMunsell(colors), 
  parseMunsell(colors, returnLAB=TRUE), 
  stringsAsFactors = FALSE
)

# make the figure
plot(B ~ A, data=x, type='n', asp=0.75, xlab='A-coordinate', ylab='B-coordinate', ylim=c(-25, 45), main='Hue Order per TN #2\nCIELAB Colorspace')
abline(h=0, v=0, lty=2)
points(B ~ A, data=x, pch=15, col=x$hex, cex=4.5)
text(x$A, x$B, labels = hue.labels, cex=0.75)

2.3 Color Contrast Charts

Inspect the data used to generate the figures with the additional argument returnData=TRUE.

Basic chart, single reference hue.

contrastChart(m = '7.5YR 4/3', hues = c('7.5YR'))

Threshold color chips by \(\Delta{E_{00}}\).

contrastChart(m = '7.5YR 4/3', hues = c('7.5YR'), thresh = 15)

Multiple reference hues.

contrastChart(m = '7.5YR 4/3', hues = c('5YR', '7.5YR', '10YR'))

Split chips by contrast class, single reference hue.

contrastChart(m = '7.5YR 4/3', hues = c('7.5YR'), style = 'CC')

Split chips by contrast class, multiple reference hues. This works better when run through latticeExtra::useOuterStrips().

fig <- contrastChart(m = '7.5YR 4/3', hues = c('5YR', '7.5YR'), style = 'CC')
useOuterStrips(fig, strip = strip.custom(bg=grey(0.85)), strip.left = strip.custom(bg=grey(0.85)) )

2.4 Example 1

# load full LUT from {aqp}
data(munsell)

# subset some common colors
hues <- c('10R', '2.5YR', '5YR', '7.5YR', '10YR', '2.5Y', '5Y')
x <- subset(munsell, subset=value %in% 3:5 & chroma %in% 1:3 & hue %in% hues)

# convert into an ordered factor, according to hue position
huePos <- huePosition(x=NULL, returnHues = TRUE)
x$hue <- factor(x$hue, levels=huePos, ordered = TRUE)
x$hue <- droplevels(x$hue)


# ok
table(x$hue)

## 
##   10R 2.5YR   5YR 7.5YR  10YR  2.5Y    5Y 
##     9     9     9     9     9     9     9

nrow(x)

## [1] 63

# convert into hex notation for plotting
x$color <- munsell2rgb(x$hue, x$value, x$chroma)

## Notice: converting hue to character

x$munsell <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)

Arrange colors in a familiar manner.

xyplot(value ~ chroma | hue, data=x,
       xlim=c(0.5, 3.5), ylim=c(2.5, 5.5), 
       scales=list(alternating=3, x=list(at=1:3), y=list(at=3:5)),
       as.table=TRUE, strip=strip.custom(bg='grey'),
       subscripts=TRUE,
       panel=function(xx, yy, subscripts, ...) {
         panel.grid(-1, -1)
         panel.points(xx, yy, col=x$color[subscripts], pch=15, cex=6)
       }
)

Evaluate \(\Delta~E_{00}\) via {farver} pacakge.

# dE00
# pair-wise distances are only present in the upper triangle
# lower triangle is 0
d <- compare_colour(x[, c('L', 'A', 'B')], from_space='lab', white_from = 'D65', method='cie2000')
# copy color codes for convenience
dimnames(d) <- list(x$munsell, x$munsell)

# convert to dist object, which expects distances in the lower triangle
d.dist <- as.dist(t(d))
# hierarchical clustering
h <- as.phylo(as.hclust(diana(d.dist)))
# nMDS
mds <- MASS::sammon(d.dist, trace = FALSE)

# remove 0's in lower triangle and diagonal
d[lower.tri(d, diag = TRUE)] <- NA
hist(d)

# double-check that ordering is correct
# yes
table(h$tip.label == x$munsell)

## 
## TRUE 
##   63

par(mar=c(0.5,0.5,1,0.5), bg='black', fg='white')

plot(h, label.offset=0.55, edge.color='white', tip.color='white', cex=0.66, direction='downwards', font=1)
tiplabels(pch=15, cex=1.5, col=x$color, offset=0.23)
title('Divisive Hierarchical Clustering', col.main='white', line=0)
mtext(text = expression("Distance metric:"~Delta*E['00']), side = 3, adj = 0, line=-0.5)

plot(h, edge.color='white', tip.color='white', cex=0.5, direction='downwards', show.tip.label=FALSE)
tiplabels(pch=15, cex=1.5, col=x$color, offset=0.25)
tiplabels(text=x$value, cex=0.55, frame='none', col='white', offset=0.23)
title('Divisive Hierarchical Clustering', col.main='white', line=0)
mtext(text = expression("Distance metric:"~Delta*E['00']), side = 3, adj = 0, line=-0.5)

par(mar=c(0.5,0.5,1.5,0.5), bg='black', fg='white')

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], x$value, cex=0.75)
mtext(text = expression("Distance metric:"~Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell value', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], paste0(x$value, '/', x$chroma), cex=0.66)
mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell value/chroma', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], x$hue, cex=0.66)
mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell hue', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

plot(mds$points, type='n', axes=FALSE)
abline(h=0, v=0, col=grey(0.85), lty=3)
points(mds$points, pch=15, col=x$color, cex=3.5)
text(mds$points[, 1],  mds$points[, 2], x$munsell, cex=0.66)
mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
mtext(text = 'Munsell color', side = 1, adj = 1, line=-1)
title('nMDS', col.main='white', line=0.75)

2.4.1 Nearly Imperceptible Differences

# locate color pairs which are at the edge of perceptable differences
idx <- which(d < 4)
row.idx <- row(d)[idx]
col.idx <- col(d)[idx]
low.dE00 <- cbind(dimnames(d)[[1]][row.idx], dimnames(d)[[2]][col.idx])

# not so helpful viz of pairs
swatchplot(list(
  A=parseMunsell(low.dE00[, 1]),
  B=parseMunsell(low.dE00[, 2])
))

# compute median delta-E00 by Munsell chip
l <- list()
for(i in seq_along(x$munsell)) {
  m <- x$munsell[i]
  d.i <- na.omit(c(d[i, ], d[, i]))
  l[[i]] <- data.frame(m=m, stat=median(d.i), stringsAsFactors = FALSE)  
}

# flatten
E00.summary <- do.call('rbind', l)

# check
kable(E00.summary[grep('10YR', E00.summary$m), ], digits=2)

	m	stat
10	10YR 3/1	11.15
11	10YR 3/2	10.64
12	10YR 3/3	11.92
13	10YR 4/1	10.45
14	10YR 4/2	10.35
15	10YR 4/3	11.37
16	10YR 5/1	11.82
17	10YR 5/2	11.58
18	10YR 5/3	13.32

# combine chip summary with source data
z <- merge(x, E00.summary, by.x='munsell', by.y='m', all.x=TRUE)

# chip size is proportional to median delta-E00
xyplot(value ~ chroma | hue, data=z,
       xlim=c(0.5, 3.5), ylim=c(2.5, 5.5), 
       scales=list(alternating=3),
       as.table=TRUE, strip=strip.custom(bg='grey'),
       subscripts=TRUE,
       panel=function(xx, yy, subscripts, ...) {
         panel.grid(-1, -1)
         panel.points(xx, yy, col=z$color[subscripts], pch=15, cex=z$stat * 0.4)
       }
)

2.5 Color Chips Arranged by Perceptual Differences

# subset some common colors
hues <- c('7.5YR')
x <- subset(munsell, subset=value %in% 3:8 & chroma %in% c(1,2,3,4,6,8) & hue %in% hues)

# convert into hex notation for plotting
x$color <- munsell2rgb(x$hue, x$value, x$chroma)
x$munsell <- sprintf("%s/%s", x$value, x$chroma)

# dE00
# pair-wise distances are only present in the upper triangle
# lower triangle is 0
d <- compare_colour(x[, c('L', 'A', 'B')], from_space='lab', white_from = 'D65', method='cie2000')
# copy color codes for convenience
dimnames(d) <- list(x$munsell, x$munsell)

# convert to dist object, which expects distances in the lower triangle
d.dist <- as.dist(t(d))

# nMDS
mds <- MASS::sammon(d.dist, trace = FALSE)

# rotate 90 degrees CCW
# to roughly follow Munsell color book page layout
# column-order
m <- matrix(
  c(0, -1, 
    1, 0), 
  byrow = FALSE, ncol = 2
)

# apply transformation
mds.trans <- mds$points %*% m

par(mar=c(2.5, 0.5, 3, 0.5), bg='black', fg='white')

plot(mds.trans, type='n', axes=FALSE, asp=1)
abline(h=seq(-25, 25, 5), v=seq(-25, 25, 5), col='white', lty=3)
points(mds.trans, pch=22, bg=x$color, cex=6)
text(mds.trans[, 1],  mds.trans[, 2], paste0(x$value, '/', x$chroma), cex=0.66)
# mtext(text = expression(Delta*E['00']), side = 3, adj = 0)
# mtext(text = '7.5YR', side = 1, adj = 1, line=-1)
title('Perceptual Distances\n7.5YR Page', col.main='white', line=0.75)
mtext(text = expression(Delta*E['00']%->%nMDS%->%90*degree~CCW~rotation), side = 1, adj = 0)

Map perceptual arrangement to original representation in color book.

p <- procrustes(mds.trans, x[, c('value', 'chroma')], scale = TRUE, symmetric = TRUE)

par(mar=c(2.5, 0, 3, 0), bg='black', fg='white')

plot(p, pch=22, kind=0, axes=FALSE, xlab='', ylab='')
points(p, display='rotated', pch=22, bg=x$color, cex=5)
points(p, display='target', pch=21, bg=x$color, cex=4)
lines(p, type='arrows', length=0.1)
text(p, display='rotated', cex=0.66, labels=paste0(x$value, '/', x$chroma))
title('Perceptual vs Color Book Arrangement', col.main='white', line=1)
title(sub='7.5YR Page', col.sub='white', line=0)
legend('bottomright', legend=c('color book', 'perceptual'), pch=c(22, 21), pt.bg=parseMunsell('7.5YR 4/6'), pt.cex=2, bty='n', inset=c(0.25, 0))

2.6 Moist vs. Dry Colors

x <- fetchOSD('musick', colorState = 'dry')
y <- fetchOSD('musick', colorState = 'moist')

m1 <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)
m2 <- sprintf("%s %s/%s", y$hue, y$value, y$chroma)

cc <- colorContrast(m1, m2)

colorContrastPlot(m1, m2, labels=c('Dry', 'Moist'), d.cex = 0.9)

Tinker with plot settings and optimize for dark background.

x <- fetchOSD('Leefield', colorState = 'dry')
y <- fetchOSD('Leefield', colorState = 'moist')

m1 <- sprintf("%s %s/%s", x$hue, x$value, x$chroma)
m2 <- sprintf("%s %s/%s", y$hue, y$value, y$chroma)

par(bg='black', fg='white')
colorContrastPlot(m1, m2, labels=c('Dry', 'Moist'), printMetrics = TRUE, d.cex = 1.25, col.cex = 1.5, label.cex = 2, label.font = 2)

2.7 Differences Between Adjacent Chips

# fresh start
data(munsell)

# for now, most commonly used hue pages, but expanded to full range in value / chroma
x <- subset(munsell, subset = hue %in% c('10R', '2.5YR', '5YR', '7.5YR', '10YR', '2.5Y', '5Y'))

# create all pair-wise combinations
cols <- paste0(x$hue, ' ', x$value, '/', x$chroma)
z <- combn(cols, 2)

# convert to 2-column data.frame
# ~ 674,541 combinations!
z <- data.frame(t(z), stringsAsFactors = FALSE)

## TODO: make this a new function

# borrowed from aqp::colorContrast
m1.pieces <- parseMunsell(z[[1]], convertColors = FALSE)
m2.pieces <- parseMunsell(z[[2]], convertColors = FALSE)

# convert to value and chroma to numeric
m1.pieces[[2]] <- as.numeric(m1.pieces[[2]])
m1.pieces[[3]] <- as.numeric(m1.pieces[[3]])
m2.pieces[[2]] <- as.numeric(m2.pieces[[2]])
m2.pieces[[3]] <- as.numeric(m2.pieces[[3]])

# difference in number of hue chips, clock-wise, as specified in:
# https://www.nrcs.usda.gov/wps/portal/nrcs/detail/soils/ref/?cid=nrcs142p2_053569
dH <- abs(huePosition(m1.pieces[[1]]) - huePosition(m2.pieces[[1]]))
# difference in number of value chips
dV <- abs(m1.pieces[[2]] - m2.pieces[[2]])
# difference in number of chroma chips
dC <- abs(m1.pieces[[3]] - m2.pieces[[3]])

## ^^^^ need a new function

# just adjacent chips
# ~ 11,498 combinations
idx <- which(dH %in% c(0,1) & dV %in% c(0,1) & dC %in% c(0,1))
z <- z[idx, ]

# compute color contrast metrics on adjacent chips
d <- colorContrast(z[[1]], z[[2]])

# summaries
quantile(d$dE00)

##        0%       25%       50%       75%      100% 
##  0.000000  4.765694  7.994458  9.797715 17.252535

# encode combinations
d$code <- sprintf("%s-%s-%s", d$dH, d$dV, d$dC)

# compute select quantiles of dE00 by combination
adj <- lapply(split(d, d$code), function(i) {
  qq <- round(quantile(i$dE00, probs = c(0.05, 0.5, 0.95)))
  data.frame(
    dH = i$dH[1],
    dV = i$dV[1],
    dC = i$dC[1],
    q05 = qq[1],
    q50 = qq[2],
    q95 = qq[3]
  )
    
}
)

adj <- do.call('rbind', adj)

# all combinations
knitr::kable(adj)

	dH	dV	dC	q05	q50	q95
0-0-1	0	0	1	0	2	4
0-1-0	0	1	0	6	8	10
0-1-1	0	1	1	7	9	11
1-0-0	1	0	0	1	4	6
1-0-1	1	0	1	2	5	7
1-1-0	1	1	0	6	10	13
1-1-1	1	1	1	6	10	13

# just changes along a single dimension 
knitr::kable(adj[row.names(adj) %in% c('0-0-1', '1-0-0', '0-1-0'), ])

	dH	dV	dC	q05	q50	q95
0-0-1	0	0	1	0	2	4
0-1-0	0	1	0	6	8	10
1-0-0	1	0	0	1	4	6

2.8 Example 2

library(rms)
library(tactile)

data(munsell)
x <- subset(munsell, subset=value %in% 3:8 & chroma %in% 1:6 & hue %in% c('10R', '2.5YR', '5YR', '7.5YR', '10YR', '2.5Y', '5Y'))

cols <- paste0(x$hue, ' ', x$value, '/', x$chroma)
z <- combn(cols, 2)

z <- data.frame(t(z), stringsAsFactors = FALSE)
d <- colorContrast(z[[1]], z[[2]])

kable(head(d), row.names = FALSE)

m1	m2	dH	dV	dC	dE00	cc
10R 3/1	10R 3/2	0	0	1	5.266520	Faint
10R 3/1	10R 3/3	0	0	2	9.026377	Distinct
10R 3/1	10R 3/4	0	0	3	11.849272	Distinct
10R 3/1	10R 3/5	0	0	4	14.044517	Prominent
10R 3/1	10R 3/6	0	0	5	15.819853	Prominent
10R 3/1	10R 4/1	0	1	0	8.811914	Faint

bwplot(cc ~ dE00, data=d, 
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Soil Color Contrast Class\nhue 10R-2.5Y, value 3-8, chroma 1-6',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

bwplot(dV ~ dE00, data=d, 
       subset = dH < 1 & dC < 1,
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Change in Munsell Value\nhue 10R-2.5Y, value 3-8, chroma 1-6',
       sub = 'Constant Hue and Chroma',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

bwplot(dH ~ dE00, data=d,
       subset = dV < 1 & dC < 1,
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Change in Munsell Hue\nhue 10R-2.5Y, value 3-8, chroma 1-6',
       sub = 'Constant Value and Chroma',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

bwplot(dC ~ dE00, data=d, 
       subset = dV < 1 & dH < 1,
       par.settings=tactile.theme(), 
       xlab=expression(Delta*E['00']),
       main='CIE Delta-E00 vs. Change in Munsell Chroma\nhue 10R-2.5Y, value 3-8, chroma 1-6',
       sub = 'Constant Hue and Value',
       scales=list(x=list(tick.number=10)), 
       panel=function(...) {
         panel.grid(h=3, v=-1)
         panel.bwplot(...)
       })

idx.F <- which(d$cc == 'Faint' & d$dE00 > 12)
idx.D <- which(d$cc == 'Distinct' & d$dE00 > 20)
idx.P <- which(d$cc == 'Prominent' & d$dE00 < 20)


kable(head(d[idx.F, ]), row.names = FALSE)

m1	m2	dH	dV	dC	dE00	cc
10R 3/1	10R 5/1	0	2	0	18.68596	Faint
10R 3/1	10R 5/2	0	2	1	18.94273	Faint
10R 3/2	10R 5/1	0	2	1	19.83064	Faint
10R 3/2	10R 5/2	0	2	0	18.72647	Faint
10R 3/2	10R 5/3	0	2	1	18.78994	Faint
10R 3/2	5Y 3/2	6	0	0	13.85876	Faint

# plot(dE00 ~ dH, data=d)
# plot(dE00 ~ dV, data=d)
# plot(dE00 ~ dC, data=d)

dd <- datadist(d)
options(datadist="dd")

# there is a little bit of curvature
(m <- ols(dE00 ~ rcs(dH, 3) + rcs(dV, 3) + rcs(dC, 3), data=d))

## Linear Regression Model
##  
##  ols(formula = dE00 ~ rcs(dH, 3) + rcs(dV, 3) + rcs(dC, 3), data = d)
##  
##                     Model Likelihood    Discrimination    
##                           Ratio Test           Indexes    
##  Obs   31626    LR chi2    100532.26    R2       0.958    
##  sigma2.3709    d.f.               6    R2 adj   0.958    
##  d.f.  31619    Pr(> chi2)    0.0000    g       12.552    
##  
##  Residuals
##  
##       Min       1Q   Median       3Q      Max 
##  -9.01452 -1.78228 -0.08724  1.68269 10.22195 
##  
##  
##            Coef   S.E.   t      Pr(>|t|)
##  Intercept 6.1933 0.0490 126.37 <0.0001 
##  dH        0.4303 0.0236  18.24 <0.0001 
##  dH'       1.1100 0.0325  34.14 <0.0001 
##  dV        4.4707 0.0260 171.90 <0.0001 
##  dV'       3.7974 0.0290 130.99 <0.0001 
##  dC        0.4921 0.0260  18.92 <0.0001 
##  dC'       0.7916 0.0290  27.31 <0.0001 
## 

# (m <- ols(dE00 ~ dH + dV + dC, data=d))

plot(Predict(m, dH=NA, dV=0, dC=0))

plot(Predict(m, dH=0, dV=NA, dC=0))

plot(Predict(m, dH=0, dV=0, dC=NA))

plot(Predict(m, dH=NA, dV=0:2, dC=0))

plot(Predict(m, dV=NA, dH=0:2, dC=0))

plot(Predict(m, dC=NA, dV=0:2, dH=0))

plot(summary(m, dH=c(0,1), dV=c(0,1), dC=c(0,1)))

plot(summary(m))

plot(anova(m), what='partial R2')

plot(nomogram(m))

For next time, how do you use these tables?

print(nomogram(m, dH=c(0, 1, 2, 3)))

## Points per unit of linear predictor: 2.535368 
## Linear predictor units per point   : 0.39442 
## 
## 
##  dH Points
##  0  0     
##  1  1     
##  2  3     
##  3  6     
## 
## 
##  dV Points
##  0    0   
##  1   12   
##  2   27   
##  3   49   
##  4   74   
##  5  100   
## 
## 
##  dC Points
##  0   0    
##  1   1    
##  2   3    
##  3   7    
##  4  11    
##  5  15

---

This document is based on aqp version 1.29.